Below is a list of useful packages for making maps in r.
library(dplyr) #for data manipulation and stuff
library(sf) #to read and process vector data
library(terra) # to read and process raster data
library(tidyterra) #useful tools for raster mapping
library(tmap) #visualization package
library(leaflet) #interactive mapping package
library(bcdata) #access to the bc data catalog
library(bcmaps) #access to common spatial data via bc data catalog
library(ggplot2) # visualization package
library(lidR) #used to read, process and visualize lidar data
library(rnaturalearth) # used to download data from natural earth
library(rnaturalearthdata) # used to download data from natural earth
library(ggspatial) #useful functions for adding scale bars and north arrows
library(geodata) #download geospatial data
library(patchwork) #create inset maps and complex figures
This is my favourite and probably the best book for all things geospatial in R: https://r.geocompx.org/
These tutorials us raster data but is a good introduction https://rspatial.org/
This is a workshop for using Lidar data. It is more advanced but useful if you would like to learn: https://tgoodbody.github.io/lidRtutorial/
https://earthengine.google.com/
https://opendata.nfis.org/mapserver/nfis-change_eng.html
library(MetBrewer)
library(viridis)
library(fishualize)
library(RColorBrewer)
library(scico)
#Example colours
MetBrewer::met.brewer(name = "VanGogh2")
fishualize::fishualize(option= "Oncorhynchus_mykiss")
scico::scico_palette_show()
You just received a grant to study the effects of elevation gradients on salmon sightings. Like any good scientist you are going to start your research by making a study area map.
First we need to download or read all of the required data
##Download data directly from URL,using BC maps, geodata, and read in files from hardrive
#Great Lakes
url <- "https://www.sciencebase.gov/catalog/file/get/530f8a0ee4b0e7e46bd300dd"
dest <- "./_data/great_lakes.zip"
if (!file.exists(dest)) {
download.file(url, destfile = dest, mode = "wb")
}
unzip_dir <- "./_data/great_lakes"
unzip(dest, exdir = unzip_dir)
great_lakes <- list.files(unzip_dir, pattern = "hydro.*\\.shp$", recursive = T, full.names = T) %>%
purrr::map(vect) %>%
purrr::map(aggregate) %>%
vect() %>%
project("epsg:3347")
data_dir = "./_data/"
# Get canadian provinces
cad <- geodata::gadm("Canada", level = 1, path = data_dir)
#get other countries
world <- geodata::world(path =data_dir) %>%
project("epsg:3347")
# read in vancouver island shapfile
study_area = vect("./_data/vancouver-island_1235.geojson")
#get BC
bc = bc_bound()%>%
vect()
Next lets start building our map. We will begin by adding all of the layers using the geom_spatvector function.
ggplot()+
geom_spatvector(data = world)+
geom_spatvector(data = cad ) +
geom_spatvector(data = great_lakes) +
geom_spatvector(data = bc )
Notice how even though we added all the layers we are still to zoomed out to see our study area? We can use the function coord_sf() to set limits in our map. But first we will programmatically determine the map extents in r
cad_ext <- cad %>%
project(world) %>%
ext()
ggplot()+
geom_spatvector(data = world)+
geom_spatvector(data = cad ) +
geom_spatvector(data = great_lakes) +
geom_spatvector(data = bc ) +
coord_sf(xlim = c(cad_ext$xmin, cad_ext$xmax),
ylim = c(cad_ext$ymin, cad_ext$ymax))
Good Job! Now we can add some color and a red box that shows the location of our study area.
# read in study area data
study_area = vect("./_data/vancouver-island_1235.geojson")
study_ext = study_area %>%
project("epsg:3347")%>%
ext()
study_box = study_ext %>%
vect("epsg:3347")
ggplot() +
geom_spatvector(data = world, fill = "#a3aeb6") +
geom_spatvector(data = cad , fill = "#c6c9ce") +
geom_spatvector(data = great_lakes, fill = "#97afb9") +
geom_spatvector(data = bc, fill = "#e5e5e5") +
geom_spatvector(data = study_box, col = "red", fill = "#00000000", linewidth = 1) +
coord_sf(xlim = c(cad_ext$xmin, cad_ext$xmax),
ylim = c(cad_ext$ymin, cad_ext$ymax))
Now lets make some final adjustments to the theme and background colors and assign the map to an object called canadamap so that we can come back to it later.
canadamap = ggplot() +
geom_spatvector(data = world, fill = "#a3aeb6") +
geom_spatvector(data = cad , fill = "#c6c9ce") +
geom_spatvector(data = great_lakes, fill = "#97afb9") +
geom_spatvector(data = bc, fill = "#e5e5e5") +
geom_spatvector(data = study_box, col = "red", fill = "#00000000", linewidth = 1) +
coord_sf(xlim = c(cad_ext$xmin, cad_ext$xmax),
ylim = c(cad_ext$ymin, cad_ext$ymax))+
theme_bw() +
theme(panel.background = element_rect(fill = "#97afb9"),
axis.ticks = element_blank(),
axis.text = element_blank(),
plot.background = element_blank())
canadamap
Congratulations you have now made a map of Canada that shows the location of a study area on Vancouver island. In the next section we will create a map of Vancouver island.
Lets start by getting a bit more data.
#download and process the digital elevation model for the study area
study_area_sf = study_area%>%
st_as_sf()
dem = cded_terra(aoi = study_area_sf)
plot(dem)
plot(st_geometry(study_area_sf), add = T)
Good you now have a digital elevation model for the study area. Notice how the elevation model extends outside of the study area. Lets fix that.
study_area_proj = study_area%>%
project(dem) #rpoject the study area shapefile to the same crs as the dem
demcrop = mask(dem,study_area_proj)
plot(demcrop)
Much better. Now that we have an elevation model we can move into ggplot to further customize the map.
ggplot()+
geom_spatraster(data = demcrop)
Lets adjust the colour and customize the legend.
ggplot()+
geom_spatraster(data = demcrop )+
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA)
Good! Now lets make the elevation really pop by creating a hillshade for a 3d effect and adjusting the theme.
slope <- terrain(demcrop, "slope", unit = "radians")
aspect <- terrain(demcrop, "aspect", unit = "radians")
hill <- shade(slope, aspect, 30, 270)
# normalize names
names(hill) <- "shades"
# Hillshading, but we need a palette
pal_greys <- hcl.colors(1000, "Grays")
ggplot()+
geom_spatraster(data = hill, alpha = 1) +
scale_fill_gradientn(colors = pal_greys,
na.value = NA,
guide = "none") +
ggnewscale::new_scale_fill()+
geom_spatraster(data = demcrop, alpha = 0.75 )+
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA)+
theme_void()
Next, because we got rid of the graticule and background using theme_void() we are going to add a north arrow, scale bar and some final layout changes.
ggplot()+
geom_spatraster(data = hill, alpha = 1) +
scale_fill_gradientn(colors = pal_greys,
na.value = NA,
guide = "none") +
ggnewscale::new_scale_fill()+
geom_spatraster(data = demcrop, alpha = 0.7 )+
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA,
breaks = c(0,1000,2000))+
theme_void()+
theme(
legend.position = "inside",
legend.position.inside = c(0.4, 0.20))+
ggspatial::annotation_north_arrow(location = "tr", which_north = "true",
pad_x = unit(1, "in"), pad_y = unit(0.5, "in"),
height = unit(2,"cm"), width = unit(2,"cm"),
style = ggspatial::north_arrow_nautical(
fill = c("grey40", "white"),
line_col = "grey"
))+
ggspatial::annotation_scale(
location = "bl",
pad_x = unit(1.7, "in"), pad_y = unit(0.3, "in"),
bar_cols = c("grey40", "white"),
text_family = "Arial",
width_hint = 0.3
)
Now lets bring all of the techniques together and create a publication ready study area map show study plot locations on Vancouver island.
We will begin by creating plot locations.
study_plots = st_sample(study_area_sf, 10)%>%
st_as_sf()%>%
mutate(Plot_label = 1:10)
#lets plot them on an interactive map
plots_latlon = study_plots %>%
st_transform(crs = 4269)
leaflet()%>%
addTiles()%>%
addMarkers(data = plots_latlon, popup = as.character(plots_latlon$Plot_label))
Lets start by adding our study area and study plots to a map
ggplot()+
geom_spatvector(data = study_area)+
geom_sf(data = study_plots, aes(col = "Study Points"))
This time lets start by making it a bit more pretty right off the bat.
ggplot()+
geom_spatvector(data = study_area)+
geom_sf(data = study_plots, aes(col = "Study Points"), size = 2)+
scale_color_manual(values = "grey20", labels = "Study Sites")+
labs(color = NULL)+
theme_bw()+
theme(legend.key=element_rect(fill=NA))
Now its starting to look ok. Lets add a few more layers and reall make it sparkle.
ggplot()+
geom_spatvector(data = bc, fill = "#e5e5e5")+
geom_spatraster(data = hill, alpha = 1) +
scale_fill_gradientn(colors = pal_greys,
na.value = NA,
guide = "none") +
ggnewscale::new_scale_fill()+
geom_spatraster(data = demcrop, alpha = 0.75) +
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA) +
geom_spatvector(data = study_area, fill = NA)+
geom_sf(data = study_plots, aes(col = "Study Points"), size = 2)+
scale_color_manual(values = "grey20", labels = "Study Sites")+
labs(color = NULL)+
theme_bw()+
theme(legend.key=element_rect(fill=NA))
Oh no, because we added new layers to the map we have changed the extent.Lets go ahead and create a study area extent object to fix this problem.
study_ext_bc = study_area %>%
project(bc)%>% #the map take the crs of the first layer
ext()
ggplot()+
geom_spatvector(data = bc, fill = "#e5e5e5")+
geom_spatraster(data = hill, alpha = 1) +
scale_fill_gradientn(colors = pal_greys,
na.value = NA,
guide = "none") +
ggnewscale::new_scale_fill()+
geom_spatraster(data = demcrop, alpha = 0.75) +
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA) +
geom_spatvector(data = study_area, fill = NA)+
geom_sf(data = study_plots, aes(col = "Study Points"), size = 2)+
scale_color_manual(values = "grey20", labels = "Study Sites")+
labs(color = NULL)+
coord_sf(xlim = c(study_ext_bc$xmin, study_ext_bc$xmax),
ylim = c(study_ext_bc$ymin, study_ext_bc$ymax))+
theme_bw()+
theme(legend.key=element_rect(fill=NA))
Good job! Now lets make some final layout changes
study_area = ggplot()+
geom_spatvector(data = bc, fill = "#e5e5e5")+
geom_spatraster(data = hill, alpha = 1) +
scale_fill_gradientn(colors = pal_greys,
na.value = NA,
guide = "none") +
ggnewscale::new_scale_fill()+
geom_spatraster(data = demcrop, alpha = 0.75) +
scale_fill_hypso_c(palette = "usgs-gswa2", name = "Elevation (m)",
guide = guide_colourbar(direction = "horizontal",
order = 2,
title.position = "top"), na.value = NA) +
geom_spatvector(data = study_area, fill = NA)+
geom_sf(data = study_plots, aes(col = "Study Points"), size = 2)+
scale_color_manual(values = "grey20", labels = "Study Sites")+
labs(color = NULL)+
coord_sf(xlim = c(study_ext_bc$xmin, study_ext_bc$xmax),
ylim = c(study_ext_bc$ymin, study_ext_bc$ymax))+
theme_bw()+
theme(
legend.position = "inside",
legend.position.inside = c(0.85, 0.85),panel.background = element_rect(fill = "#97afb9"),
legend.background = element_rect(fill = NA),
legend.text = element_text(hjust = 0.5, size = 9),
legend.title = element_text(hjust = 0.5),
legend.key=element_rect(fill=NA),
legend.spacing.y = unit(0,"cm"),
)
study_area
Now look at all of that blank space in the corner of the map. If only we had a map of canada laying around so that people from all over the world reading our article could figure out where this is….. Oh wait…
#we can use the patchworks package to create an inset map
Final = study_area + inset_element(canadamap, 0,0,0.4,0.4)
Final
Lastly, lets export the map with correct sizing. This can be journal specific but I generally default to the elsevier sizing. For this map we want it to be a full page width and roughly a 3rd of a page tall.
ggsave("./StudyAreaMap.jpeg", plot = Final, dpi = 300, units = "mm", width = 185, height = 150)